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Abstract 

We study numerically the applicability of the effective-viscosity approach for sim- 
ulating the effect of gravitational instability (GI) in disks of young stellar objects 
with different disk-to-star mass ratios ^. We adopt two a-parameterizations for the 
effective viscosity based on Lin & Pringle (1990) and Kratter et al. (2008) and com- 
pare the resultant disk structure, disk and stellar masses, and mass accretion rates 
with those obtained directly from numerical simulations of self-gravitating disks 
around low-mass (M^, 1.0 Mq) protostars. We find that the effective viscosity 
can, in principle, simulate the effect of GI in stellar systems with ^ ^ 0.2 — 0.3, thus 
corroborating a similar conclusion by Lodato & Rice (2004) that was based on a dif- 
ferent a-parameterization. In particular, the Kratter et al's a-parameterization has 
proven superior to that of Lin & Pringle's, because the success of the latter depends 
crucially on the proper choice of the a-parameter. However, the a-parameterization 
generally fails in stellar systems with ^ ^ 0.3, particularly in the Class and Class 

I phases of stellar evolution, yielding too small stellar masses and too large disk-to- 
star mass ratios. In addition, the time-averaged mass accretion rates onto the star 
are underestimated in the early disk evolution and greatly overestimated in the late 
evolution. The failure of the a-parameterization in the case of large ^ is caused by a 
growing strength of low-order spiral modes in massive disks. Only in the late Class 

II phase, when the magnitude of spiral modes diminishes and the mode-to-mode 
interaction ensues, may the effective viscosity be used to simulate the effect of GI 
in stellar systems with ^ 0.3. A simple modification of the effective viscosity that 
takes into account disk fragmentation can somewhat improve the performance of a- 
models in the case of large ^ and even approximately reproduce the mass accretion 
burst phenomenon, the latter being a signature of the early gravitationally unsta- 
ble stage of stellar evolution (Vorobyov & Basu, 2006). However, further numerical 
experiments are needed to explore this issue. 
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1 Introduction 

It has now become evident that circumstellar disks are prone to the devel- 
opment of gravitational instability in the early stage of stellar evolution (e.g. 
Laughlin & Bodenheimer, 1994; Vorobyov & Basu, 2005; Vorobyov & Basu, 
2006; Krumholz et al., 2007; Kratter et al., 2008; Stamatellos & Whitworth, 
2008). The non-axisymmetric spiral structure resulting from GI produces grav- 
itational torques, which are negative in the inner disk and positive in the outer 
disk and help limit disk masses via radial transport of mass and angular mo- 
mentum (Laughlin et al., 1997; Adams et al., 1989; Vorobyov & Basu, 2005). 
Even in the late cvohition. weak gravitational torques associated with low- 
amplitude density perturbations in the disk can drive mass accretion rates 
typical for T Tauri stars (Vorobyov & Basu, 2008). 

The fact that gravitational torques trigger mass and angular momentum redis- 
tribution in the disk makes them conceptually similar to viscous torques, which 
are believed to operate in a variety of astrophysical disks. An anticipated ques- 
tion is then whether the Gl-induced transport in circumstellar disks can be 
imitated by some means of effective viscosity? The answer depends on whether 
mass and angular momentum transport induced by gravitational torques is of 
global or local nature and this issue still remains an open question. 

Lin & Pringle 1987; 1990 were among the first to suggest that the transport 
induced by GI could be described within a viscous framework and parameter- 
ized the effect of GI using the following formulation for the effective kinematic 
viscosity 



where Q is the angular speed, Cg is the sound speed, Qt = CsO/(7rGE) is 
the Toomre parameter, E is the gas surface density, and is the critical 
value of Qt at which the disk becomes unstable against nonaxisymmetric 
instability. The dimensionless number «lp represents the efficiency of mass 
and angular momentum transport by GI. It is evident that Lin & Pringle's 
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(1990) parameterization is actually that of Shakura & Sunyacv's (1973), with 
a-parameter modified by the term {Ql/Q\ — 1) to describe the effect of GI. 

Laughlin & Rozyczka (1996) have compared the evolution of a thin, self- 
gravitating protostellar disk using two-dimensional hydrodynamic simulations 

with the evolution of a one-dimensional, axisymmetric viscous disk using the 
model of Lynden-Bell & Pringle (1974). They suggest the following effective 
a-parameter as a modification to the usual Shakura & Sunyaev a-model 



(2) 



where the surface density E and angular velocity Vt are in specific units defined 
in Laughlin & Rozyczka (1996), and C, and j3 are defined as 

C-0.325 - 2.25 i \ , 

/5 = 0.875 -4.347, 

where M* and are the star and disk masses, respectively, and 7 is the poly- 
tropic constant. By considering models with different masses of protostellar 
disks and different values of Qtj they concluded that while their parameteriza- 
tion works better than that of Shakura & Sunyaev, the a-prescription fails in 
relatively massive disks with ^ 0.5 M* due to the presence of global modes 
which are not tightly wound. Laughlin & Rozyczka's (1996) parameterization 
also performs badly in low-mass disks with Afj 0.2 M*, overestimating the 
actual radial mass transport due to gravitational torques. In the intermediate 
mass regime, however, their parameterization can reproduce the actual gas 
surface density profiles with a good accuracy. The appearance of dimensional 
units in the dimensionless coefficient cklr restricts the applicability of equa- 
tion (2) - one has to cither use the Laughlin & Rozyczka's (1996) system of 
units or redefine equation (2) according to the adopted units. 

Lodato & Rice (2004) have studied the apphcability of the viscous treatment 
of the circumstellar disk evolution. They have modelled the evolution of self- 
gravitating circumstellar disks with masses ranging from 5% to 25% of the 
star using adiabatic equation of state with 7 = 5/3 and cooling that removes 
energy on timescales of 7.5 orbital periods. By calculating the gravitational 
and Reynolds stress tensors and assuming that the sum of these stresses is 
proportional to the gas pressure, they estimated the effective a-parameter 
associated with gravitational and Reynolds stresses to be «g^R < 0.06 in their 
numerical simulations. By further assuming that their disks quickly settle in 
thermal equilibrium (when the rate of viscous energy dissipation is balanced 
by radiative cooling) , they proposed the following expression for the saturated 
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value of the a-parameter (see also Gammie, 2001) 



eq 



(ilnr 



1 



7(7 - l)tcool^' 



(5) 



where tcooi = 7.50~^ is a characteristic disk cooling time. For the parameters 
adopted in Lodato & Rice's (2004) numerical simulations, = 0.05. A good 
agreement between numerically derived Q;g3 and cteq allowed them to conclude 
that the viscous treatment of self-gravity is justified in disks with disk-to-star 
mass ratios ^ < 0.25 and, perhaps, even in more massive disks (Lodato & 
Rice, 2005). 

Recently, Kratter et al. (2008) have suggested a modification to the usual 
Shakura & Sunyaev a-model based on the previous numerical simulations of 
Laughlin & Rozyczka (1996), Lodato & Rice 2004; 2005, and Gammie (2001). 
Their a-parameter invoked to simulate the effect of GI (ctci) consists of two 
components: the "short" component ccshort and "long" component cciong 

«short + «longj > (6) 



where 



"short = max 



along = max 



0.14 



1^ 



-l)(l-/.)^■^^o 



1.4 X 10-3(2 - Qt) 



,0 



(7) 
(8) 



where // is the disk-to-total mass ratio. The "short" component aghort differs 
from the effective a-parameter suggested by Lin & Pringle's equation (1) only 
in a mild dependence. The "short" and "long" components are meant to 
represent the different wavelength regimes of the gravitational instability ex- 
pected to dominate at different values of \i and Qt- We note that Kratter et 
al. also assumes Qt = max((5T, 1)- 

The use of the a-model as a proxy for Gl-induccd transport has been chal- 
lenged by Balbus & Papaloizou (1999), who argue that the energy flux of 
self-gravitating disks is not reducible to a superposition of local quantities 
such as the radial drift velocity and stress tensor - the essence of an a-disk. 
Instead, extra terms of non-local nature are present that mostly invalidate the 
a-approach. Their conclusion is corroborated by recent numerical simulations 
of collapsing massive cloud cores by Krumholz et al. (2007), who flnd that 
gravitational instability in embedded circumstellar disks with masses of order 
50% that of the central star is dominated by the m — \ spiral mode induced by 
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the SLING instability (Adams et al., 1989; Shu et al, 1990). This mechanism 
is global and enables mass and angular momentum transport on dynamical 
rather than viscous timescales. 

Circumstellar disks form in different physical environments and go through 
different phases of evolution so that it is quite likely that there is no unique 
answer on whether Gl-induccd transport can be described by some means 
of effective viscosity. Indeed, the early embedded phases of disk evolution 
(Class O/I) are substantially influenced by an infalling envelope, both through 
the mass deposition (Vorobyov & Basu, 2006; Kratter et al., 2008) and enve- 
lope irradiation (Matzner & Levin, 2005; Cai et al., 2008). Younger Class O/I 
disks are usually more massive than in the older Class II ones (Vorobyov, 
2009). As a result, gravitational instabilities in Class O/I disk may be domi- 
nated by low-order (m < 2), large-amplitude global modes (see e.g. Krumholz 
et al., 2007), which are likely to invalidate the viscous approach. On the other 
hand. Class II disks evolve in relative isolation and settle into a steady state 
that is characterized by low-amplitude, high-order modes m > 3 (e.g. Lodato 
& Rice, 2004; Vorobyov & Basu, 2007). Such modes tend to produce more 
fluctuations and cancellation in the net gravitational torque on large scales, 
thus making the viscous approach feasible. 

The mass and angular momentum transport in self-gravitating disks is diffi- 
cult to deal with analytically due to a kaleidoscope of competing spiral modes. 
Furthermore, long-term multidimensional numerical simulations of the evolu- 
tion of circumstellar disks involving an accurate calculation of disk self-gravity 
are usually very computationally intensive. On the other hand, the theory of 
viscous accretion disks is fairly well developed (e.g. Pringlc, 1981) and is rel- 
atively easy to deal with numerically. This motivated many authors to use 
the viscous approach to mimic the effect of gravitational instability when 
studying the long-term evolution of circumstellar disks (e.g. Lin & Pringle, 
1990; Nakamoto & Nakagawa, 1995; Hueso & Guillot, 2005; DuUemond et al., 
2006; Kratter et al., 2008). It is therefore important to know if and when the 
viscous approach is applicable. Some work in this direction has already been 
done by Lodato & Rice 2004; 2005 and Cossins et al. (2009), who considered 
a short-term evolution of isolated disks with different disk-to-star mass ratios 
meant to represent different stellar evolution phases. However, they have not 
considered a long-term disk evolution due to an enormous CPU time demand 
of fully three-dimensional simulations. In the present paper we explore the 
applicability of the a-parameterization of gravitational instability along the 
entire stellar evolution sequence, starting from a deeply embedded Class 
phase and ending with a late Class II phase (T Tauri phase). Although the 
T Tauri phase of stellar evolution is likely to harbour only marginally grav- 
itationally unstable disks with associated torques of low intensity (Vorobyov 
& Basu, 2007), the disk structure may bear the imprints of the early, GI- 
dominated phase. That is why it is important to capture the main stages of 
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disk evolution altogether in one numerical simulation. We focus on the a- 
parameterizations of Lin & Pringle (1990) and Krattcr ct al. (2008) and defer 
a study of the Lodato & Rice's (2004) parameterization to a follow-up paper. 
Contrary to other studies, we run our numerical simulations of circumstellar 
disks first with self-gravity calculated accurately by solving the Poisson inte- 
gral and then with self-gravity imitated by effective viscosity. We then perform 
a detailed analysis of the resultant circumstellar disk structure, masses, and 
mass accretion rates in the both approaches. 



2 Description of numerical approach 



2.1 Main equations 



We seek to capture the main evolution phases of a circumstellar disk alto- 
gether, starting from the deeply embedded Class phase and ending with 
the T Tauri phase. This can be accomplished only by adopting the so-called 
thin-disk approximation. In this approximation, the basic equations for mass 
and momentum transport are written as (Vorobyov & Basu, 2006; Vorobyov 
& Basu, 2009) 

— = -Vp • (Ev,) , (9) 
= -V,P + E9i + E9l+{V- n)^ , (10) 



where E is the mass surface density, V = J^2: ^(^^ the vertically integrated 
form of the gas pressure P, Z is the radially and azimuthally varying vertical 
scale height, Vp = VrV + v^(f) is the velocity in the disk plane, Vp = rd/dr + 
(j)r~^d/d(j) is the gradient along the planar coordinates of the disk, and g'^'^ = 

gf'^r + g'^^(t) is the gravitational acceleration in the disk plane. The latter 
consists (in general) of two parts: that due to the disk self-gravity {g^) and 
that due to the gravity of the central star (gip . The gravitational acceleration 
gfp is found by solving for the Poisson integral (Vorobyov & Basu, 2006). The 
viscous stress tensor 11 is expressed as 



n = 2E I/, 



eff 



(v^;-^(V-^;)e), (11) 



where Vf is a symmetrized velocity gradient tensor, e is the unit tensor, and 
i/eff is the effective kinematic viscosity. The components of (V ■ 11)^ in polar 
coordinates (r, 0) can be found in Vorobyov & Basu (2009). We emphasize that 
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we do not take any simplifying assumptions about the form of the viscous stress 
tensor, apart from those imposed by the adopted thin-disc approximation. It 
can be shown (Lodato, 2008) that equation (10) can be reduced to the usual 
equation for the conservation of angular momentum of a radial annulus in the 
axisymmetric viscous accretion disc (Pringle, 1981). 

Equations (9) and (10) are closed with a barotropic equation that makes 
a smooth transition from isothermal to adiabatic evolution at E = Ecr ~ 
36.2 g cm-^: 

7' = c^E + c^E„(^)\ (12) 

where = 0.188 km s~^ is the sound speed in the beginning of numerical 
simulations and 7 = 1.4. Assuming a local vertical hydrostatic equilibrium in 
the disk, Ecr = 36.2 g cm~^ becomes equivalent to the critical number density 
10"^^ cm-3 (Masunaga & Inutsuka, 2005). 

The thin-disk approximation is an excellent means to calculate the evolution 
for many orbital periods and many model parameters. It is well justified as 
long as the aspect ratio A = A{r) of the disk vertical scale height Z to radius 
r does not considerably exceed 0.1. The aspect ratio A{r) for a Keplerian disk 
is usually approximated by the following expression 

A<CQ^M^{r)/M,, (13) 

where Md(r) is the disk mass contained within radius r and C is a constant, 
the actual value of which depends on the gas surface density distribution E 
in the disk. For a disk of constant surface density, C is equal unity. However, 
circumstellar disks are characterized by surface density profiles declining with 
radius. For the scahng E oc r'^-^ typical for our disks, C — 1/4. Adopting 
further Qt = 2.0 and M^{r)/M^ — 0.5, which are the upper limits in our nu- 
merical simulations, we obtain the maximum value of A{r) = 0.25. This anal- 
ysis demonstrates that the thin-disk approximation may be only marginally 
valid in the outer regions of a circumstellar disk, but is certainly justified in 
its inner regions where Md(r)/M^ is small. A typical distribution of the aspect 
ratio A in a disk around one solar mass star was shown in figure 7 in Vorobyov 
(2009). 

2.2 Initial conditions 

We start our numerical integration in the pre-stellar phase, which is charac- 
terized by a collapsing (flat) starless cloud core, and continue into the late 



7 



accretion phase, which is characterised by a protostar-disk system. This en- 
sures a self-consistent formation of a circumsteUar disk, which occupies the 
innermost portion of the computational grid, while the infalling envelope (in 
the embedded stage of stellar evolution) occupies the rest of the grid. 



We consider model cloud cores with mass Md = 1.5 M0, initial temperature 
T = 10 K, mean molecular weight 2.33, and the outer radius rout = 0.07 pc. 
The initial radial surface density and angular velocity profiles are characteristic 
of a collapsing axisymmetric magnetically supercritical core (Basu, 1997): 



J.2 _|_ J.2 



(14) 



r 




(15) 



where VLq is the central angular velocity. The scale length tq = kc1/{GIlo), 
where k = \/2/7i and Sq = 0.12 g cm~^. These initial profiles are characterized 
by the important dimensionless free parameter r] = ^0^0/^5 and have the 
property that the asymptotic (r ^ Tq) ratio of centrifugal to gravitational 
acceleration has magnitude V2r] (Basu, 1997). The centrifugal radius of a mass 
shell initially located at radius r is estimated to be rcf — j^/{Gm) — y/2rjr, 
where j = Qr'^ is the specific angular momentum. 

The strength of gravitational instability is expected to depend on the disk 
mass. According to Lodato & Rice (2004,2005), a viscous parameterization of 
GI is allowed in systems with the disk-to-star mass ratio ^ <^ 0.25 but may 
fail in systems with relatively more massive disks. It is therefore important to 
consider systems with different disk-to-star mass ratios. This can be achieved 
by varying the initial rate of rotation of a model cloud core, but keeping all 
other cloud core characteristics fixed. Indeed, an increase in Qq would result 
in a larger value of rj and rcf. In turn, an increase in the centrifugal radius r^f 
would result in more mass landing onto the disk rather than being accreted 
directly onto the central star (plus some inner circumsteUar region at r < 5 AU 
which is unresolved in our numerical simulations), thus raising the resultant 
disk-to-mass ratio. According to Caselli et al. (2002), velocity gradients in 
dense molecular cloud cores range between 0.5 and 6.0 km s~^ pc^^. Therefore, 
we choose three typical values for the central angular velocity of our model 
cloud cores: Q — Qo,i — 0.82 km s~^ pc"^, O, — fio,2 = 1-5 km s~^ pc~^, and 
= ^0,3 = 2.5 km s~^ pc~^. It is convenient to parameterize our models in 
terms of the ratio of rotational to gravitational energy f3, which is very similar 
in magnitude to the parameter i] introduced above. Caselli et al. (2002) report 
f3 ranging between 10~^ and 0.07 in their sample of dense molecular cloud 
cores. Our model cloud cores have Pi — 8.83 x 10~^, P2 — 2.3 x 10~^, and 
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(3^ = 8.2 X 10 ^. For the sake of conciseness, we present the results of the /3i 
and models, referring to the intermediate (32 model only where necessary. 

Equations (9), (10), (12) are solved in polar coordinates (r, cj)) on a numerical 
grid with 128 x 128 cells. We use the method of finite differences with a time- 
explicit, operator-split solution procedure. Advection is performed using the 
second-order van Leer scheme. The radial points are logarithmically spaced. 
The innermost grid point is located at Tin = 5 AU, and the size of the first ad- 
jacent cell is 0.3 AU. We introduce a "sink cell" at r < 5 AU, which represents 
the central star plus some circumstellar disk material, and impose a free in- 
fiow inner boundary condition. The sink cell is dynamically inactive but serves 
as the source of gravity, thereby influencing the gas dynamics in the active 
computational grid. The outer boundary is reflecting. The gravity of a thin 
disk is computed by directly summing the input from each computational cell 
to the total gravitational potential. The convolution theorem is used to speed 
up the summation. A small amount of artiflcial viscosity is added to the code 
to smear out shocks over one computational zone according to the usual von 
Neumann & Richtmyer's (1950) prescription, though the associated artificial 
viscosity torques were shown to be negligible in comparison with gravitational 
torques (Vorobyov & Basu, 2007). A more detailed explanation of numerical 
methods and relevant tests can be found in Vorobyov & Basu 2006; 2007. 

2.3 Three numerical models 

We consider three numerical models that are distinct in the way the right-hand 
side of equation (10) is handled. In the first model (hereafter, self-gravitating 
model or SG model), viscosity is set to zero throughout the entire evolution 
(forth term) and the system evolves exclusively via gravity of the disk and 
central star (second and third terms, respectively), as well as pressure forces 
(first term). In the second model (hereafter, Lin and Pringle model or LP 
model) , the disk self-gravity is switched off after the disk formation (no second 
term) and the subsequent evolution is governed by pressure forces (first term) , 
gravity of the central star (third term), viscosity (forth term). In particular, 
the kinematic viscosity is computed using the following representation 



where h is the disk scale height. The latter is calculated assuming the vertical 
hydrostatic equilibrium in the gravitational field of the disk and the central 
star (see Vorobyov, 2009). For numerical reasons, i/efr is set to zero if Qt drops 
accidentally below some low value, which is set in our simulations to 0.3. The 




(16) 
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third model (hereafter, KMK model) differs from the second model only in the 
way the effective viscosity is defined. More specifically, we use the following 
expression for the effective kinematic viscosity 

i^eflf = OlGlCsh, (17) 

where ckgi is determined from equation (6). Following Kratter et al. (2008), 
we also set Qt = max(QTi !)• In Section 6, we consider a modification to the 
KMK model that attempt to deal with the fragmentation regime at Qt < 
1.0. All the three numerical models start from identical initial conditions as 
described in Section 2.2 



3 Cloud cores with low rates of rotation. 

In this section we consider model cloud cores that are described by the ratio of 
rotational to gravitational energy f3 = f3i = 0.88 x 10~^. This value is chosen 
to represent dense cloud cores with low rates of rotation, as inferred from the 
measurements of Caselli et al. (2002). 

3.1 The LP model 

Equation (16) indicates that the LP model has two free parameters: cklp and 
Qcr- It is therefore our main purpose to determine the values of q;lp, us- 
ing which the LP model reproduces best the exact solution provided by the 
SG model. We have found that the LP model depends only weakly on Qcr, 
as soon as its value is kept near 1.5. In this study we set Qcr — 1-7. We vary 
ctLp in wide limits, starting from lO^'^ and ending with 0.5. Figures 1-3 show 
the disk radial structure obtained in the SG model (solid lines) and LP model 
(dashed lines). More specifically, the top panels show the radial gas surface 
density distributions, while the middle and bottom panels show the radial pro- 
files of the gas temperature and Toomre parameter Q^, respectively. It may 
not be entirely consistent to calculate Qt in the non-self-gravitating model. 
Nevertheless, this quantity may serve as an indicator of how far the stability 
properties of the LP model deviate from the exact solution. Three character- 
istic times based on the age of the central star in the SG model have been 
chosen: t = 0.2 Myr (left row), t = 1.0 Myr (middle row), and t = 2.0 Myr 
(right row). The left row represents an evolution stage when about 25% of 
the initial cloud core material is still contained in the infalling envelope, while 
the middle and right rows represent a stage when the envelope has almost 
vanishes. Three values of the a-parameter have been selected for the presen- 
tation: q;lp = 10-^ (Fig. 1), ctLP = 0.01 (Fig. 2), and cclp = 0.25 (Fig. 3). 
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Fig. 1. Gas surface density (top row), gas temperature (middle row), and 
Q-parameter (bottom row) as a function of radius in the cklp = 10"^ model (dashed 
lines) and SG model (solid lines). Three distinct evolutionary times since the forma- 
tion of the central protostar are indicated in each column. The dotted lines show the 
MMSN gas surface density profile (top row) and a gas temperature profile typical 
for circumstellar disk (middle row) as inferred by AW05. 

The dotted lines show the minimum mass solar nebular (MMSN) density pro- 
file S = 1.7 X I0'^(r/1AU)~^'^ (Hayashi, 1985) and a disk radial temperature 
profile T = 180(r/lAU)-°-^ inferred by Andrews & Williams (2005) (hereafter 
AW05) from a large sample of YSO in the Taurus- Aurigae star-forming region. 

It is instructive to review the main disk properties obtained by the SG model. 
An accurate determination of disk masses in numerical simulations of collaps- 
ing cloud cores is not a trivial task. Self-consistently formed circumstellar disks 
have a wide range of masses and sizes, which are not known a priori. In addi- 
tion, they often experience radial pulsations in the early evolution phase, which 
makes it difficult to use such tracers as rotational support against gravity of 
the central star. However, numerical and observational studies of circumstellar 
disks indicate that the disk surface density is a declining function of radius. 
Therefore, we distinguish between disks and inf ailing envelopes using a critical 
surface density for the disk-to-envelope transition, for which we choose a fidu- 
cial value of Str = 0.1 g cm~^. This choice is dictated by the fact that densest 
starless cores have surface densities only slightly lower than the adopted value 
of Etr- In addition, our numerical simulations indicate that self-gravitating 
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Fig. 2. The same as Fig. 1 only for the cklp = 10 ^ model. 

disks have sharp outer edges and the gas densities ol order 0.01 — 0.1 g cm~^ 
characterize a typical transition region between the disk and envelope. 

As the solid lines in Figs 1-3 indicate, most of the self-gravitating disk is 
characterized by a power-law surface density distribution declining with ra- 
dius as r~^'^. This slope is also predicted by the MMSN hypothesis (Hayashi, 
1985). However, our obtained gas surface densities are almost a factor of 10 
greater than those of the MMSN throughout the entire disk evolution. This 
feature is an important property of self-gravitating disks (see also Vorobyov 
& Basu, 2009). It makes easier giant planets to form, because planet forma- 
tion scenarios seem to require gas densities at least a few times larger than 
those of the minimum-mass disk (Ida & Lin, 2004; Boss, 2001). The radial 
gas temperature profiles indicate that the self-gravitating disk is somewhat 
colder than the typical disk in Taurus- Aurigae region (AW05). This is par- 
ticularly true in the late evolution, though large deviations from the typical 
profile toward colder disk are also present in the AW05 sample. We point out 
that irradiation by a central source can raise the disk temperature and se- 
riously affect the disk propensity to fragmentation in the inner disk regions 
(e.g. Matzner & Levin, 2005). Unfortunately, this effect is difficult to take into 
account self-consistently in polytropic disks due to the lack of detailed treat- 
ment of cooling and heating. Hotter disk can be obtained in our numerical 
simulations by choosing a larger value for the ratio of specific heats 7 in the 
barotropic equation of state (12). We discuss hotter disks in Section 5. We 
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Fig. 3. The same as Fig. 1 only for the q;lp = 0.25 model, 
also note that our self-gravitating disk exhibits a near-Keplerian rotation. 

Finally, the time behaviour of the Toomre parameter Qt in a self-gravitating 
disk warrants some attention. It is evident that the self-gravitating disk reg- 
ulates itself near the boundary of gravitational stabiUty, with values of Qt 
lying in the 1.7 — 2.0 range. This important property breaks down when a 
sufficient amount of physical viscosity (i.e. due to turbulence) is present in 
the disk (see e.g. Vorobyov & Basu, 2009). We note that the Q-parameter 
in the early evolution of a self-gravitating disk may occasionally drop below 
1.7 — 2.0 in some parts of the disk. These episodes are usually associated with 
disk fragmentation and formation of dense clumps, which have masses of up to 
10-20 Jupiters, sizes of several AU, number densities of up to 10^^ — 10^^ cm~^ 
and are pressure supported against their own gravity. Most of these clumps 
are quickly driven onto the central protostar by gravitational torques from 
spiral arms. This phenomenon causes a burst of mass accretion (Vorobyov & 
Basu, 2005; Vorobyov & Basu, 2006) and is very transient in nature (the burst 
itself takes less than 100 yr). A small number of the clumps may be flung into 
the outer regions where they disperse, most likely due to a combined action 
of tidal forces, differential rotation, and insufficient numerical resolution (our 
grid is logarithmically spaced in the radial direction). 

We now proceed with comparing the disk structure in the SG and LP models. 
Figures 1-3 reveal that both the q;lp = 10~^ model and cklp = 0.25 model 
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fail to accurately reproduce the radial structure of the disk obtained in the 
SG model. More specifically, the cilp = 10~^ model yields too dense and hot 
disks throughout the entire evolution. The disagreement is particularly strong 
in the early disk evolution {t — 0.2 Myr), when the radial gas surface density 
distribution becomes substantially shallower than that predicted by the exact 
solution (oc r~^'^). The obtained disk radii in the q;lp = 10~^ model are smaller 
by a factor of 2 — 3 than those found in the SG model. The Toomre parameter 
is also considerably smaller than that of the self-gravitating disk. In summary, 
the LP model with cklp as small as 10~^ fails to provide an acceptable fit to 
the exact solution. 

When we turn to large values of the a-parameter q;lp = 0.25 (Fig 3), the 
corresponding LP model seems to yield density and temperature profiles that 
are in acceptable agreement with those of the exact solution only in the very 
early phase of disk evolution {t — 0.2 Myr). Even in this early stage, the disk 
size is severely underestimated. Furthermore, the late evolution sees a strong 
mismatch between the disk structure in the LP and SG models. We conclude 
that large values of q;lp ^ 0.25 may be marginally acceptable in the early, 
strongly gravitationally unstable phase of disk evolution, but cannot be used 
to simulate the effect of self-gravity on long time scales of order of several 
Myr. 

The ttLP = 10~^ model (dashed lines in Fig. 2) appears to provide the best fit 
to the exact solution. Although in the early disk evolution (t — 0.2 Myr) the 
simulated gas surface densities and temperatures are somewhat larger than 
those provided by the SG model, in the late evolution (t 1 Myr) they 
are almost indistinguishable from the exact solution throughout most of the 
disk. We conclude that the LP model can be rather successful in reproducing 
the effect of gravitational instability in star-disk systems formed form slowly 
rotating cloud cores, provided that the value of cclp lies close to 10~^. 

It is interesting to compare disk and stellar masses obtained in the LP models 
with those of the SG model. The left column in Fig. 4 shows the stellar mass 
(top), disk mass (middle), and disk-to-star mass ratio ^ (bottom) for the SG 
and LP models. More specifically, the solid, dashed, dash-dotted, and dotted 
lines present data for the SG model, Olp = 10~^ model, Olp = 10~^ model, 
and q;lp = 0.25 model, respectively. We use Sent = 0.1 g cm^^ for the disk-to- 
envelope transition. The horizontal axis shows time elapsed since the formation 
of the central star. The disk forms at t — 0.14 Myr ^ , when the central object 
has accreted approximately 60% of the initial cloud core mass Md = 1.5 Mq. 
The disk-to-star mass ratio in the SG model never exceeds ^ = 17%. It is 



In fact, the disk forms earlier but its evolution is unresolved in the inner 5 AU due 
to the use of a sink cell in our numerical simulations. However, the mass contained 
in this inner 5 AU is negligible compared to the rest of the resolved disk. 
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Fig. 4. Stellar masses (top row), disk masses (middle row), and disk-to-star mass 
ratios (bottom row) in the LP model (left column) and KMK model (right column). 
In particular, the dashed, dash-dotted, and dotted lines in the left column present 
data for the aLP = 10"^ model, olp = 10^^ model, and olp = 0.25 model, respec- 
tively, while the dashed line in the right column shows data for the KMK model. 
The solid lines in both columns correspond to the SG model. 



evident that the ctLp = 10~^ model yields the masses and disk-to-star mass 
ratios that agree best with those derived in the SG model. The high- viscosity 
LP model (cklp = 0.25) predicts too large values for the stellar mass and too 
low values for the disk mass, whereas the low-viscosity LP model (cklp = 10""') 
does it vice versa. This example nicely illustrates the sensitivity of the LP 
model to the choice of olp- 

Finally, we calculate the instantaneous mass accretion rates through the inner 
disk boundary as M = —2Tir--aJlv^^ where r^^ = 5 AU, and fj. is the radial 
gas velocity at r^^. There is evidence that accretion on to the star is a highly 
variable process in the early embedded stage of disk evolution. Observations 
indicate that, in addition to young stellar objects (YSOs) with mass accretion 
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Fig. 5. Time-averaged mass accretion rates in the SG model (solid line), q;lp = 10~^ 
LP model (dashed line), olp = 10~^ LP model (dash-dotted line), and qlp = 0.25 
LP model (dotted line). The insert shows the instantaneous mass accretion rate 
versus time in the SG model. 



rates similar to those predicted by Shu (1977), there is a substantial populace 
of YSOs with the sub-Shu accretion rates M < 10^^ Mq yr^^ (Enoch et al., 
2009) and a small number of super-Shu accretors with M > 10~^ Mq yr~^. 
Furthermore, numerical simulations show that the disk in the early embedded 
phase becomes periodically destabilized due to the mass deposition from an 
infalling envelope (Vorobyov & Basu, 2005; Vorobyov & Basu, 2006; Boley, 
2009) . The resulted gravitational torques drive excess mass in the form of 
dense clumps on to the central star, thus producing the so-called burst phe- 
nomenon (Vorobyov & Basu, 2005; Vorobyov & Basu, 2006). This phenomenon 
can only be captured approximately by any model that mimics mass transport 
in self-gravitating disks via an effective viscosity (see e.g. Kratter et al., 2008). 
In Section 6, we try to reproduce the burst behaviour within the framework 
of the KMK model by modifying the definition of the a-parameter. 

In order to smooth out the bursts and facilitate a comparison between the 
SG and LP models, we calculate the mean mass accretion rates (M) by time- 
averaging the instantaneous rates over 2 x lO'' yr. Figure 5 shows (M) versus 
time for the SG model (solid line), cclp = 10~^ LP model (dashed line), cclp = 
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10~^ LP model (dash-dotted line), and olp = 0.25 LP model (dotted line). 
The horizontal axis shows time elapsed since the formation of the central star. 
The burst phenomenon is illustrated in the insert to Fig. 5, which shows the 
instantaneous mass accretion rates in the SG model as a function of time. It 
is seen that the time-averaged mass accretion rates are identical before the 
disk formation {t < 0.14 Myr) but become distinct soon afterward. Even after 
time-averaging, some substantial variations in the mass accretion rate of the 
SG model are still visible. It is evident that the q;lp = 10""^ LP model greatly 
underestimates (M) in the early 0.4 Myr, while considerably overestimating it 
in the subsequent evolution. The high- viscosity cklp = 0.25 LP model also seem 
to produce larger accretion rates than those of the exact solution, especially 
in the late phase. The «lp = 10^^ LP model, again, demonstrates better 
agreement with the exact solution, though somewhat underestimating (M) 
after 1 Myr. 



3.2 The KMK model 



In this section we investigate the efficiency of the Kratter et al.'s (2008) a- 
parameterization in simulating the effect of GI in circumstellar disks. The 
KMK model described by equations (6)-(8) is more useful and flexible than 
the LP model because the former has no free parameters^ and it includes 
a mild dependence on the disk-to-total mass ratio /x. This ratio is expected 
to depend on the initial rate of rotation of a molecular cloud core and may 
considerably vary along the sequence of stellar evolution phases, and therefore 
could be an important ingredient for successful modeling of the effect of Gl 
in circumstellar disks. For a more detailed discussion we refer the reader to 
Kratter et al. (2008). 

Figure 6 shows radial profiles of the gas surface density (top), temperature 
(middle), and Toomre parameter (bottom) obtained in the KMK model (dashed 
lines) and in the SG model (solid lines). The horizontal axis, as usual, shows 
time elapsed since the formation of the central star. Three characteristic stel- 
lar ages, as indicated in each row, are chosen for the presentation. It is evident 
that the KMK model shows a satisfactory performance, particularly in the 
late evolution stage. In the early evolution, the predicted disk surface den- 
sities and temperatures are slightly larger than those of the SG model, but 
the difference is not significant. Comparing Figs 2 and 6 one can see that the 
KMK model reproduces the exact solution to the same extent and accuracy as 
the best q;lp = 10~^ LP model. However, the supremacy of the KMK model 
is obvious - it has no free parameters to adjust. We believe that the KMK 



^ In fact, the value of the critical Toomre parameter is set to 1.3 in the Kratter 
et al.'s approach. 
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Fig. 6. Gas surface density (top row), gas temperature (middle row), and 
Q-parameter (bottom row) as a function of radius in the KMK model (dashed lines) 
and SG model (solid lines). Three distinct evolutionary times since the formation 
of the central protostar arc indicated in each column. The dotted lines show the 
MMSN gas surface density profile (top row) and a typical gas temperature profile 
(middle row), as inferred fir circumstellar disks by AW05 



model owes its success to the use of a two-component a-parameter, aci- 

The right column in Fig. 4 shows the stellar mass (top), disk mass (middle), 
and disk-to-star mass ratio ^ obtained using the KMK model (dashed lines) 
and SG model (solid lines). This figure demonstrates that the KMK model 
yields the disk and stellar masses that arc in good agreement with those of 
the SG model. We conclude that both the LP and KMK models can, in prin- 
ciple, reproduce the radial structure of sclf-gravitating disks, time-averaged 
accretion rates, and masses in star/disk systems formed from slowly rotating 
cloud cores with /3 ^ 1.0 x 10^^. Such systems are characterized by disk-to-star 
mass ratios not exceeding ^ = 0.2. The a-parameterization in such systems 
appears to be justified. The agreement with the exact solution is somewhat 
modest in the early evolution but improves considerably in the late evolution 
as the strength of gravitational instability declines. In the case of LP model, 
the a-parameter has to be set to cklp ~ 0.01. 
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Fig. 7. Prom top to bottom: gas surface density, gas temperature, angular velocity, 

and Q-parameter as a function of radius in the KMK model (dashed lines) and 
SG model (solid lines) for the (3-^ = 8.2 x 10~^ case. Three distinct evolutionary 
times since the formation of the central protostar are indicated in each column. 
The dotted lines show the MMSN gas surface density profile (first row) , typical gas 
temperature profile (second row) as inferred in circumstellar disks by AW05, and a 
Kepler rotation law (third row). 

4 Cloud cores with intermediate and high rates of rotation. 



In this section we consider cloud cores that are characterized by the ratios of 
rotational to gravitational energy P2 — 2.3 x 10~^ and Ps — 8.2 x 10~^, with 
most attention being concentrated on the latter case. Since the KMK model 
has proven superior in comparison to the LP model, we focus on the former 
one, summarizing the main results for the LP model where necessary. Cloud 
cores with high rates of rotation are expected to yield disks with large disk- 
to-star mass ratios. Our motivation is then to determine the extent to which 
the KMK model can reproduce the exact solution in systems with disk-to-star 
mass ratios considerably larger than that of the previous section, ^ ~ 0.17. 

Figure 7 compares the disk radial structure obtained in the SG model (solid 
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lines) and KMK model (dashed lines) for the (3^ case. The layout of the figure 
is exactly the same as in Fig. 1, but we have also plotted the radial distribution 
of the angular velocity Q in the third row from top. The dotted line in this 
row shows a Kepler rotation law, fl oc r~^-^. 

Both the approximate and exact solutions are characterized by similar slopes 
of the gas surface density S oc r^^'^, gas temperature T oc r^°'^, and angular 
velocity fl oc r~^^^. One can see that the KMK model reproduces the radial 
gas surface densities and temperatures but yields too low angular velocities 
and the Toomre parameter, especially in the early evolution. This implies that 
the KMK model underestimates the mass of the central star. 

The failure of the KMK model to accurately reproduce stellar masses is illus- 
trated in Fig. 8, which show the stellar masses (top row), disk masses (middle 
row), and disk-to-star mass ratios (bottom row) for the SG model (solid hnes) 
and KMK model (dashed lines). In particular, the left column corresponds 
to the (32 case, while the right column presents the jS^ case. The horizontal 
axis shows time elapsed since the formation of the central star. To compare 
masses along the sequence of stellar evolution phases, we need an evolutionary 
indicator to distinguish between Class 0, Class I, and Class II phases. We use 
a classification of Andre et al. (1993), who suggest that the transition between 
Class and Class 1 objects occurs when about 50% of the initial cloud core 
is accreted onto the protostar-disk system. The Class II phase is consequently 
defined by the time when the infalling envelope clears and its total mass drops 
below 10% of the initial cloud core mass M^. The vertical dotted hnes mark 
the onset of Class I (left line) and Class II (right line) phases. A general trend 
of the KMK model to underestimate the stellar masses and to overestimate 
the disk-to-star mass ratios is clearly seen. 

To quantify this mismatch between the models, we calculate time-averaged 
stellar masses (M*), disk masses (Md), and disk-to-star mass ratios in 
each major stellar evolution phase. The resulted values are listed in Table 1. 
It is evident that the mismatch between the SG model and KMK model is 
particularly strong in the Class and Class I phases. For instance, the /^a 
KMK model yields (^) = 1.31 in the Class I phase, which is almost a factor of 3 
larger than the corresponding value for the SG model. The disagreement in the 
Class II phase is in general less intense than in the Class and Class I phases. 
Somewhat surprisingly, disk masses in the KMK model differ insignificantly 
from those of the SG model. The q;lp = 10~^ model shows a very similar 
behaviour. 

Our numerical simulations demonstrate that the extent to which the KMK 
model departs from the SG model is direct proportional to the value of (3, and 
hence to the disk-to-star mass ratio. This is not surprising. Disks in systems 
with greater ^ are more gravitationally unstable, and the viscous approach is 
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Fig. 8. Time evolution of stellar masses (top row), disk masses (middle row), 
and disk-to-star mass ratios (bottom row) in the KMK model (dashed lines) and 
SG model (solid lines). The left and right columns correspond to the P2 = 2.3 x 10~^ 
and = 8.2 x 10~^ cases, respectively. The horizontal axis shows time since the 
formation of the central star. 

Table 1 

SG model versus KMK model for P2 = 2.3 x IQ-'^ and P3 = 8.2 x lO"'^ 



phase (M,)^, {Md)/3, 



Class 


0.39/0.33 


0.09/0.12 


0.23/0.33 


0.29/0.17 


0.14/0.17 


0.49/0.97 


Class I 


0.81/0.56 


0.28/0.33 


0.34/0.58 


0.77/0.31 


0.42/0.42 


0.54/1.33 


Class II 


1.05/0.98 


0.42/0.44 


0.4/0.46 


0.91/0.71 


0.55/0.64 


0.6/0.94 



All mean masses are in M©. The slash differentiates between the SG model (left) 
and KMK model (right). 
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expected to fail in massive disks, which are dominated by global spiral modes 
of low order (see e.g. Lodato & Rice, 2005; Krumholz et al., 2007). We illustrate 
this phenomenon by calculating the global Fourier amplitudes (GFA) defined 
as 
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where is the disc mass and raise is the disc's physical outer radius. The in- 
stantaneous GFA show considerable fluctuations and we have to time-average 
them over 2 x 10^ yr in order to produce a smooth output. The time evolution 
of the time-averaged GFA (log units) is shown in Fig. 9 for the /3i = 0.88 x 10~^ 
SG model (top), /32 = 2.3 x IQ-^ SG model (middle), and = 8.2 x IQ-^ SG 
model (bottom). Each color type corresponds to a mode of specific order, as 
indicated in the legend 

The time behaviour of the GFA is indicative of two qualitatively different 
stages in the disk evolution. In the early stage {t ^ 1.0 Myr), a clear segre- 
gation between the modes is evident - the lower order mode dominates its 
immediate higher order neighbour in all models. In particular, the m = 1 
mode is almost always the strongest one^ . The modes also show a clear ten- 
dency to decrease in magnitude with time. In the late stage, however, this 
clear picture breaks into a kaleidoscope of low-magnitude modes competing 
for dominance with each other. These mode fluctuations are not a numerical 
noise but are rather caused by ongoing low-amplitude non-axisymmetric den- 
sity perturbations sustained by swing amplification at the disk's sharp outer 
edge. As was shown by Vorobyov & Basu (2007), self-gravity of the disk is 
essential for these density perturbations to persist into the late disk evolution. 
The density perturbations quickly disappear if self-gravity is switched off. 

The visual analysis of GFA in the early stage reveals that models with greater 
P are characterized by spiral modes of greater magnitude. For instance, the 
magnitude of the m—1 mode a,tt — 0.5 Myr in the Pi model is Ci = —1.9 dex, 
while P2 and models have Ci = —1.7 dex and Ci = —1.5 dex, respectively. 
What is more important is that the relative strength of the low-order modes 
(m < 2) is greater in models with greater /? (and greater These facts, when 
taken altogether, account for the failure of the viscous approach in systems 
with ^ ^ 0.2. The preponderance of strong low-order modes, which are global 
in nature, largely invalidates the local-in-nature viscous approach. When we 

^ In the present paper, we have ignored a possible wobbling of the central star, 
which may increase the strength of the odd modes, especially that of the m = 1 
mode (Adams et al., 1989; Shu et al., 1990). Numerical hydrodynamics simulations 
with the indirect potential in the momentum equation (10) are needed to accurately 
assess the strength of this effect. 
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Fig. 9. Global Fourier amplitudes (averaged over 2 x 10^ yr) in the /3i = 0.88 x 10~^ 
SG model (top), (^2 = 2.3 x lO'^ SG model (middle), and /?3 = 8.2 x lO'^ SG model 
(bottom). The horizontal axis shows time elapsed since the formation of the central 
star. Each line type corresponds to a mode of specific order, as indicated in the 
legend. 
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turn to the late stage at t 1 Myr, we see that all modes saturate at a 
considerably lower value of order C = —3.5 dex. Moreover, characteristic 
fluctuations tend to produce more cancellation in the net gravitational torque 
on large scales, thus making the Gl-induced transport localwise. This explains 
why the KMK model performs better in the late stage. 

Finally, we present in Figure 10 the time-averaged mass accretion rates (M) 
as a function of time for the P2 = 2.3 x 10^"^ case (top-left panel) and /^s ~ 
8.2 X 10~^ case (top- right panel). In particular, the solid and dashed lines 
correspond to the SG model and KMK model, respectively. The time-averaged 
values are obtained from the instantaneous mass accretion rates by applying 
a running average method over 2 x 10"^ yr. For comparison, the bottom row 
presents (M) versus time for the (3i = 0.88 x 10^"^ case. More specifically, 
the dashed line in the bottom-left panel corresponds to the best cklp — 10~^ 
LP model, while the same line type in the bottom-right panel presents the 
KMK model. In both bottom panels, the solid line shows the data for the SG 
model. The horizontal axis shows time elapsed since the beginning of numerical 
simulations. 

The central star forms in all models at around 0.05 Myr, when (M) increases 
sharply to a maximum value of ~ 10^^ Mq jr~^. A period of near constant 
accretion then ensues, when the matter is accreted directly from the infalling 
envelope onto the star. This stage may be very short (or even evanescent) in 
systems with high rates of rotation (top panels) due to almost instantaneous 
onset of a disk formation stage. In this stage, the matter is first accreted onto 
the disk and through the disk onto the star, and the mass accretion rate is 
highly variable. It is evident that the (32 and /^s KMK models (top row) under- 
estimate (M) in the early disk evolution, while considerably overestimating 
(M) in the late evolution. This explains why viscous models with large ^ tend 
to greatly underestimate stellar masses in the Class and Class I phases. It 
appears that viscous torques in massive disks fail to keep up with gravita- 
tional torques in the early, strongly gravitationally unstable phase but drive 
too high rates of accretion in the late, gravitationally quiescent phase. On the 
other hand, the bottom row in Fig. 10 indicates that viscous models with low 
rates of rotation and small disk-to-star mass r atios (/9 ^ 1 X 10-3 and ^ ^ 0.2) 
show a better fit to the SG model, departing from the exact solution only by 
a factor of several. 



5 Higher disk temperature 

Observations and numerical simulations suggest that circumstellar disks are 
characterized by a variety of physical conditions, including vastly different disk 
sizes, disk-to-star mass ratios, and temperature profiles (see e.g. Andrews & 
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Fig. 10. Mass accretion rates (M) (averaged over 2 x 10^ yr) versus time. The top 
row shows (M) for (52 = 2.3 x IQ-^ (top-left panel) and /^s = 8.2 x 10"^ (top-right 
panel). In particular, the solid and dashed lines correspond to the SG and KMK 
models, respectively. The bottom row presents {Ad) for the /?i = 0.88 x 10~^ case. 
More specifically, the dashed lines in the left and right panels show the data for the 
CKLP = 10~^ LP model and KMK model, respectively, while the solid lines in both 
low panels present the data for the SG model. 

Williams, 2005; Vorobyov, 2009). We have considered the effect of different 
disk-to-star mass ratios in the previous sections. However, it is also important 
to consider disks with different temperatures, since this is one of the factors 
that control the disk propensity to gravitational instability and fragmentation. 
In our polytropic approach, we can increase the disk temperature by raising the 
ratio of specific heats 7. In a real disk, the rise in 7 does not necessarily result 
in higher temperatures due to a strong dependence on cooling and heating 
terms. These effect are not considered in the present study. 

Figure 11 presents the disk radial structure for the 7 = 5/3 and = 8.2 x 10^"^ 
case. The dashed and solid lines show the data obtained in the KMK model 
and SG model, respectively. The layout of the figure is identical to that of 
Fig. 7. The comparison of Fig. 7 and Fig. 11 reveals that the 7 = 5/3 disk is 
characterized by roughly a factor of 2 greater gas temperature than that of the 
7 = 1.4 disk. In fact, the 7 = 5/3 disk has temperatures that systematically 
exceed typical temperatures inferred by AW05 for a sample of circumstellar 
disks (dotted hue in the second row). Nevertheless, the performance of the 



25 




KMK model is largely unaffected by this change in the disk temperature. As 
in the 7 = 1.4 case, the KMK model yields too low angular velocities in the 
early evolution (third row), thus underestimating stellar masses. In general, 
the temporal behaviour of disk masses, stellar masses, and disk-to-star mass 
ratios is very similar to that shown in Fig. 8 for the 7 = 1.4 disk. This is 
not unexpected. Disk and stellar masses are largely determined by the mass 
accretion rate onto the star. As Vorobyov & Basu (2009) have demonstrated, 
a factor of 2 increase in the disk temperature makes little effect on the time- 
averaged mass accretion rates, though it may alter the temporal behaviour of 
the instantaneous mass accretion rates. 



6 Disk fragmentation 



When the local Toomre paramter Qt drops below some critical value Qcr,i, 
which is often equal or close to 1.0, fragmentation in a circumstellar disk 
ensues. To account for this qualitatively different regime of disk evolution. 
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Kratter et al. (2008) have defined a critical surface density 

CsO, 



(19) 



ttGQ, 



and assumed that fragmentation depletes the disk surface density when Qt < 



According to Kratter et al. (2008), this rate is fast enough to ensure that 
Qt never dips appreciably below (5cr,f- The actual dynamics of the fragments 
is not followed, instead they are allowed to accrete on to the star at a rate 
M^j = 0.05 Mf a 

This approach, feasible in model simulations akin to those presented by Krat- 
ter ct al. (2008), is difficult to implement in numerical hydrodynamic simu- 
lations. Taking away matter and instantaneously transporting it over a large 
distance in the numerical grid may lead to numerical instabilities of unknown 
consequences. One possible way around this problem is to actually create frag- 
ments and follow their dynamics. However, when disk self-gravity is absent, 
this approach is also misleading because the dynamics of such fragments is 
largely governed by the gravitational interaction with the disk. 

On the other hand, our viscous models demonstrate that the fragmentation 
regime is sometimes achieved in the early disk evolution (see e.g. Figure 7) 
and some changes to the standard effective-viscosity approach are necessary. 
In the present work, we modify the a-paramctcr so as to mimic a possible 
increase in the efficiency of mass transport due to fragmentation. In particular, 
in the Kratter et al.'s (2008) formulation of oqi we have made the following 
modification to the short component 



if the fragmentation regime with Qt < Qcr,i is set in the disk. In the usual 
regime of Qt > QcT,f, the short component is not modified, i.e. aghort.f = ctshort- 
Wc set Qcr,f = 1-0 and n = 10, thus introducing a strong non-linearity in 
the expression for the a-parameter. We note that in the definition of ttshort 
given by Equation (7), Qt is not allowed to drop below unity by setting 
Qt — max((5T, 1-0). In practice, this means that (Qcr.f/Qx)" is the only term 
that is sensitive to fragmentation. 

Figure 12 presents a comparison of the KMK models with and without frag- 
mentation with the SG model for the case oi P — 8.2 x 10~^ and 7 = 1.4. 



Qcr,f a-t a rate 



Ef = — (E — Ecr,f)il. 



(20) 




(21) 
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Fig. 12. Radial gas surface density profiles (top row) and Toomre Q-parameters 
(middle row) in the SG model (solid lines), unmodified KMK model (dashed lines), 
and modified- for- fragmentation KMK model (dash-dot-dotted lines). The bottom 
row presents the time evolution of the stellar mass (left panel), disk mass (middle 
panel), and disk-to-star mass ratio (right panel) in all three models. 

The top and middle rows show the radial profile snapshots of the gas surface 
density and Q-parameter at three different times after the formation of the 
central star. The bottom row shows the time evolution of the stellar mass (left 
panel), disk mass (middle panel), and disk-to-star mass ratio (right panel). 
As usual, the solid and dashed lines represent the SG model and unmodified 
KMK model, respectively. The dash-dot-dotted fine shows the KMK model 
modified for fragmentation. The vertical dotted lines mark the onset of the 
Class I and Class II phases in the SG model. 

There are several interesting features in Figure 12 that deserve attention. 
When fragmentation is taken into account in the KMK model, the Q-parameter 
never dips appreciably below the critical value for fragmentation Qcr,f = 1-0, 
contrary to what was sometimes seen in the unmodified KMK model. This 
illustrates a self-regulating nature of the modification we have made to the 
standard KMK model. Furthermore, the modified KMK model reproduces 
better the exact solution, though the mismatch is still substantial. For in- 
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Fig. 13. Time evolution of instantaneous mass accretion rates in the unmodified 
KMK model (top panel), modified-for-fragmentation KMK model (middle panel), 
and SG model (bottom panel). 
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stance, the disk-to-star mass ratio in the modified KMK model is much closer 
to that of the SG model. In particular, the disk in the modified KMK model 
is almost always less massive than the star, which is in agreement with the 
expectations of the exact SG model (Vorobyov, 2009). 

Perhaps, the most important feature of the modified KMK model is a non- 
monotonic, step-like increase in the mass of the central star with time, akin 
to that of the SG model. The disk-to-star mass ratios also demonstrates tem- 
poral variations similar to those of the SG model. We emphasize that the 
unmodified KMK model lacks such behaviour. These short-term variations 
are signatures of the mass accretion burst phenomenon. To illustrate this, 
we plot the instantaneous mass accretion rates M versus time in Figure 13 
for the unmodified KMK model (top), modified KMK model (middle), and 
SG model (bottom) . The mass accretion rate in the unmodified KMK model 
exhibits only small- amplitude flickering, there is no trace of the several-orders- 
of-magnitude variability that is typical for the burst phenomenon. Much to 
our own surprise, the modified KMK model was found to have short-term 
variations in M that are similar in amplitude to those of the SG model. Al- 
though the number of such bursts in the modified KMK model is still smaller 
than in the SG model, the mere fact that, after a simple modification, the 
KMK model can reproduce the burst phenomenon and accretion variability 
by several orders of magnitude is fascinating and encouraging. 



7 Conclusions 

We have performed numerical hydrodynamic simulations of the self-consistent 
formation and long-term (2 Myr) evolution of circumstellar disks with the 
purpose to determine the applicability of the viscous a-parameterization of 
gravitational instability in self-gravitating disks. In total, we have considered 
three numerical models: the LP model that uses the Lin & Pringle's (1990) 
a-parameterization, the KMK model that adopts the Krattcr ct al.'s (2008) 
a-parameterization, and the SG model that employs no effective viscosity but 
solves for the gravitational potential directly (the exact solution). We then 
perform a detailed analysis of the resultant circumstellar radial disk structure, 
masses, and mass accretion rates in the three models for systems with different 
disk-to-star mass ratios ^. We find the following. 

(1) The agreement between the viscous a-models and the SG model depends 
on the value of ^ and deteriorates along the sequence of increasing disk- 
to-star mass ratios. In principle, the viscous a-models can provide an 
acceptable fit to the SG model for stellar systems with ^ <^ 0.2 — 0.3, 
which is in agreement with ^ < 0.25 previously reported by Lodato & 
Rice (2004) based on a different a-parameterization. 



30 



(2) The success of the Lin & Pringle's a-parameterization in systems with 
^ ^ 0.2 — 0.3 depends crucially on the proper choice of Olp. The ctLp = 
10~^ model yields an acceptable fit to the exact solution but completely 
fails for ^ 10~^ and q;lp ^ 0.25. In particular, the q;lp = 10~^ model 
yields a disk being too small in size, having too large gas surface density 
and disk-to-star mass ratio. On the other hand, the olp = 0.25 model 
drives too large accretion rates, resulting in a disk being heavily depleted 
in mass and having too small gas-to-star mass ratios already after 1.0 Myr 
of evolution. 

(3) The performance of the KMK model is comparable to or even better than 
that of the best cklp = 10~^ LP model. In addition, the former is superior 
because it has no explicit dependence on the a-parameter, and it includes 
some dependence on the disk-to-star mass ratio. 

(4) The viscous a-models generally fail in stellar systems with { ^ 0.3. They 
yield too small stellar masses and too large disk-to-star mass ratios, espe- 
cially in the early Class and Class I evolution phases (see Table 1). For 
instance, the KMK model may overestimate the disk-to-star mass ratio 
by a factor of 3 as compared to that of the SG model. The same lack of 
agreement is also seen in the time-averaged mass accretion rates (M). In 
particular, the KMK model underestimates (M) in the Class and Class 

I phases, while greatly overestimating (M) in the Class II phase. 

(5) The failure of the KMK model (and LP model) in the case of large ^ 
is related to the growing strength of low-order spiral modes in massive 
self-gravitating disks. 

(6) A simple modification to the a-parameter that takes into account disk 
fragmentation can somewhat improve the performance of the KMK model 
and even reproduce to some extent the mass accretion burst phenomenon 
(Vorobyov & Basu, 2005; Vorobyov & Basu, 2006), demonstrating the im- 
portance of the proper treatment for disk fragmentation. More numerical 
study is needed to explore this issue. 

(7) Both viscous a-models perform better in the late disk evolution (Class 

II phase) than in the early disk evolution (Class and Class I phases), 
irrespective of the value of ^. This is due to a gradual decline in the 
magnitude of the spiral modes with time, as well as due to growing mode- 
to-mode interaction. The latter tends to produce more cancellation in 
the net gravitational torque on large scales, thus making the effect of 
gravitational torques similar to that of local viscous torques. 

(8) A factor of 2 increase in the disk temperature does not noticeably affect 
the efficiency of the viscous a-models in simulating the effect of GI. 
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